Building global hexagonal grids
Petr Keil
October 26, 2016
1 The necessary libraries and stuff
library(dggridR)
library(rgdal)
library(ggplot2)
# Shapefile of global administrative areas:
admin <- readOGR(dsn="ADMIN", layer="TM_WORLD_BORDERS")## OGR data source with driver: ESRI Shapefile
## Source: "ADMIN", layer: "TM_WORLD_BORDERS"
## with 246 features
## It has 11 fields2 Create the hexagonal grids
Construction the dggs objects
# for cell area of 629,710.6441 km^2
dggs4 <- dgconstruct(type = "ISEA3H", res = 4)
# for cell area of 209,903.54803 km^2
dggs5 <- dgconstruct(type = "ISEA3H", res = 5)
# for cell area of 69,967.84934
dggs6 <- dgconstruct(type = "ISEA3H", res = 6) Create the global grids, unfortunately in the .klm format. The resulting objects are just addresess of where the .klm files area saved.
G4 <- dgearthgrid(dggs4, frame = TRUE, wrapcells = TRUE, savegrid = TRUE)
G5 <- dgearthgrid(dggs5, frame = TRUE, wrapcells = TRUE, savegrid = TRUE)
G6 <- dgearthgrid(dggs6, frame = TRUE, wrapcells = TRUE, savegrid = TRUE)3 Switching to QGIS
Here is a critical step: I was unable to read the .klm files neither with rgdal, neither with dggridR’s dg_process_kml. Maybe it will be possible on Windows, I have no idea.
This is what happens:
dg_process_kml(G4, frame=TRUE, wrapcells=TRUE)## Error in eval(expr, envir, enclos): could not find function "dg_process_kml"What I did is: I loaded the .klm files to QGIS, and saved them as esri shapefiles.
4 Reading the ESRI shapefiles to R
G4.poly <- readOGR(dsn="GRID", layer="GRIDres4")## OGR data source with driver: ESRI Shapefile
## Source: "GRID", layer: "GRIDres4"
## with 812 features
## It has 11 fields## Warning in readOGR(dsn = "GRID", layer = "GRIDres4"): Z-dimension discarded G5.poly <- readOGR(dsn="GRID", layer="GRIDres5")## OGR data source with driver: ESRI Shapefile
## Source: "GRID", layer: "GRIDres5"
## with 2432 features
## It has 11 fields## Warning in readOGR(dsn = "GRID", layer = "GRIDres5"): Z-dimension discarded G6.poly <- readOGR(dsn="GRID", layer="GRIDres6")## OGR data source with driver: ESRI Shapefile
## Source: "GRID", layer: "GRIDres6"
## with 7292 features
## It has 11 fields## Warning in readOGR(dsn = "GRID", layer = "GRIDres6"): Z-dimension discarded5 Simple plots
par(mai=c(0.1, 0.1, 0.1, 0.1))
plot(admin, col="grey", border=NA)
plot(G4.poly, add=TRUE)
plot(admin, col="grey", border=NA)
plot(G5.poly, add=TRUE)
plot(admin, col="grey", border=NA)
plot(G6.poly, add=TRUE)
6 Fancy ggplots
6.1 Function preparing the shapefile for ggplot2
This funciton thakes the hexagonal grid, and converts it into a data frame that is usable by ggplot2. It also sorts out the ugly wrapping of the hexagons around the edges.
The function is adopted from the original dg_process_kml in package dggridR. It takes shapefiles instead of the .klm files
Arguments:
map - SpatialPolygonsDataFrame or SpatialPolygons class representing the grid.
clean.hex.map <- function(map)
{
long <- NULL
group <-NULL
map@data$timestamp <- NULL
map@data$begin <- NULL
map@data$end <- NULL
map@data$altitudeMode <- NULL
map@data$extrude <- NULL
map@data$visibility <- NULL
map@data$drawOrder <- NULL
map@data$icon <- NULL
map@data$description <- NULL
map@data$tessellate <- NULL
map@data$id <- rownames(map@data)
map.points <- fortify(map, region="id")
map.df <- merge(map.points, map@data, by="id")
# Find dangerous polygons based on how many degrees of longitude they span
groups_to_wrap <- map.df %>% group_by(group) %>%
summarise(diff=max(long)-min(long)) %>%
filter(diff>180) %>% select(group)
# Adjust them so they appear on the eastern side of the map
map.df <- map.df %>%
mutate(long=ifelse(group %in% groups_to_wrap$group,
ifelse(long<0,long+360,long), long))
# Arrange polygon points so they are ordered appropriately, otherwise the results
# will not be nice, closed cells, but weird triangular thingies
clean.map <- map.df %>% arrange(group,order)
return(clean.map)
}6.2 Processing the shapefiles to data.frames for ggplot2
gg.map4 <- clean.hex.map(G4.poly)
gg.map5 <- clean.hex.map(G5.poly)
gg.map6 <- clean.hex.map(G6.poly)
gg.admin <- fortify(admin)## Regions defined for each Polygons6.3 The plots
ggplot(gg.map4, aes(y=lat, x=long, group=group)) +
theme_bw() +
geom_polygon(data=gg.admin, aes(y=lat, x=long, group=group), fill="grey") +
geom_polygon( colour="black", fill=NA)
ggplot(gg.map5, aes(y=lat, x=long, group=group)) +
theme_bw() +
geom_polygon(data=gg.admin, aes(y=lat, x=long, group=group), fill="grey") +
geom_polygon( colour="black", fill=NA)
ggplot(gg.map6, aes(y=lat, x=long, group=group)) +
theme_bw() +
geom_polygon(data=gg.admin, aes(y=lat, x=long, group=group), fill="grey") +
geom_polygon( colour="black", fill=NA)
7 The end
sessionInfo()## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dggridR_0.1.11 dplyr_0.5.0 ggplot2_2.1.0 rgdal_1.1-9
## [5] sp_1.2-3
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.7 knitr_1.13 magrittr_1.5 munsell_0.4.3
## [5] colorspace_1.2-6 lattice_0.20-34 R6_2.1.2 stringr_1.0.0
## [9] plyr_1.8.3 tools_3.3.1 grid_3.3.1 gtable_0.2.0
## [13] DBI_0.4-1 htmltools_0.3.5 lazyeval_0.2.0 yaml_2.1.13
## [17] digest_0.6.9 assertthat_0.1 tibble_1.1 formatR_1.4
## [21] evaluate_0.9 rmarkdown_0.9.6 labeling_0.3 stringi_1.0-1
## [25] scales_0.4.0